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Abstract 

We consider the diffusion of a non-relativistic heavy quark of fixed mass M, in a 
one-dimensionally expanding and strongly coupled plasma using the AdS / CFT duality. 
The Green's function constructed around a static string embedded in a background 
with a moving horizon, is identified with the noise correlation function in a Langevin 
approach. The (electric) noise decorrelation is of order 1/T(r) while the velocity de- 
correlation is of order MD(t)/T(t). For MD > 1, the diffusion regime is segregated 
and the energy loss is Langevin-like. The time dependent diffusion constant D{t) 
asymptotes its adiabatic limit 2/-k\/~\T(t) when t/tq = (1/37/oTo) 3 where 770 is the 
drag coefficient at the initial proper time tq. 



1 Introduction 



The quark-gluon plasma (QGP) created in Relativistic Heavy Ion Collisions at RHIC is 
believed to be strongly coupled [1]. The AdS/CFT correspondence [2] has proven to be a 
useful tool for addressing issues of a strongly coupled plasma albeit in the limit of a large 
number of colors N c and strong gauge coupling A = g\ M N c . A number of nonperturbative 
properties in linear response have been recently addressed via the gravity dual calculation in 
AdS 5 black hole geometry [3]. The transport results bear some relevance to the QGP plasma 
at RHIC. Much less is known about the time-dependent evolution of the strongly coupled 
plasma. There have been suggestions that the fireball in Relativistic Heavy Ion Collision 
(RHIC) can be explained from a dual gravity point of view [4, 5]. 

In order to model the expanding plasma in a gravity set up, the use of a moving black hole 
was suggested in [4]. In [6, 7] (hereon refer to as JP), it was shown that the moving horizon 
black hole geometry can be generated by assuming an asymptotycally expanding perfect fluid 
in a boost invariant setting. Using conformal invariance and energy-momentum conservation 
together with holographic renormalization [8], JP constructed the bulk geometry from the 
perfect fluid boundary data. This metric was extended to the case of shear viscosity [11, 12], 
R-charge [13] and to a 3 dimensional non-isotropic setting [14]. An exact background with 
isotropic expansion was worked out in [15]. 

In this paper we study the diffusion of a non-relativistic heavy quarks in one dimensionally 
expanding plasma using the JP metric. Heavy quark diffusion have been studied in a static 
black hole background in various ways [16, 17, 18, 19, 20, 21, 22]. Here, we follow the 
suggestion in [20] in the static case, by analyzing the momentum fluctuations of a heavy 
quark in an expanding plasma, and use it to estimate the drag force and diffusion rate. We 
will use a generalized Langevin equation to assess the diffusion rate. Unlike [21] we suggest 
that the green function of string fluctuation should be identified with the correlator of the 
fluctuation force rather than the total force in a Langevin approach. 

The basic object of this procedure is the gauge-invariant electric-electric force decorrela- 
tion, which we will calculate using the AdS/CFT duality. The field theory dual (operator) 
to a quark displacement (£(i)|„=o) is the (gauge invariant) colored force acting on a heavy 
quark [20, 21, 22]. According to the AdS/CFT correspondence the generating function in 
field theory should be related with the classical action by 




(1) 
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whose second derivative gives us the symmetrized Wightman function 

G(h,t 2 ) = I(F(t!)F(t 2 ) + F{t 2 )F{t 1 )) . 



(2) 



It is related to the retarded Green's function [24] 

G(u>) = -coth^ImG^) , (3) 

in momentum space. We will calculate Gr(uj) from the Nambu-Goto action of the fluctu- 
ating string in the JP metric following the way proposed in [20] (section 2). We note that 
this decorrelator applies both to heavy and light fundamental quarks, therefore the decorre- 
lation time is a measure of how the gluon rescattering against external fundamental probes 
decorrelate. 

In section 2 we start with the JP metric in [6] and derive equations of motion of the 
transverse string fluctuation of a heavy quark. We also calculate the retarded Green's func- 
tion, which is translated to an electric force-force correlator (or more precisely decorrelator) 
in section 3. In section 4, we compute the momentum fluctuation(broadening) of a heavy 
quark. By arguing that a generalized Langevin equation applies for the case,we deduce the 
heavy quark diffusion property in the expanding and cooling plasma. Our conclusions are 
in section 5. 



2 String fluctuation in Janik-Peschanski metric 

We consider a heavy particle moving in an expanding plasma. Let tp = 1/T be the colored 
electric force decorrelation time (see below) and let to — MD/T be the diffusion time of a 
particle with diffusion constant D and mass M. For times tp < t < tp, o r MD > 1, the force 
decorrelation can be treated as instantaneous. Thus, the velocity decorrelation of a heavy 
quark in a strongly coupled plasma can be followed generically by a Langevin description. 
For an expanding plasma the description involves a generalized Langevin of the form 

f t+V (t)p(t)=F(t), (4) 

with both the drag coefficient rj(t) and the diffusion constant D(t) time dependent. The 
idea here is that the force-force decorrelator can be calculated from the microscopic physics 
and the drag coefficient i](t) is related to the force- force decorrelator by a non-equilibrium 
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relation [27, 28]. 

For the drag coefficient, we need to calculate the correlators in a frame which moves 
with the particle. So it is enough to consider a particle moving with the expanding plasma. 
For a one dimensional expansion, such a comoving frame was introduced by Hwa [9] and by 
Bjorken [10] in terms of the rapidity y and the proper time r of the comoving coordinate 
system, r and y are related to the time coordinate x°, and logitudinal coordinate x 3 , through 
x° = t cosh y and x 3 = r sinh y . The Minkowski metric in terms of r, y can be written as 

ds 2 = -dr 2 + r 2 dy 2 + dx 2 ± , (5) 

where x± = {x 1 ,^ 2 } are the transverse coordinates. 

The gravity dual of the Hwa- Bjorken system was worked out in [6, 7] and the metric can 
be written in the following suggestive form 



R 2 

ds 2 = tL 

z A 



(l + O 



4^2 

dr 2 + (1 + v 4 ) {r 2 dy 2 + dx 2 ± ) + dz 2 



(6) 



where z is the fifth coordinate in AdSs and v is a scaling variable defined by 

v = ~ 7~~TT , e = i(7rT ) 4 , (7) 
(t/t ) 3 4 

Since the horizon is located at v = 1 or z ~ r 1//3 , the black hole horizon may be considered 
as moving away from the boundary. We remark that the solution is valid only for late times 
(asymptotic solution). 

The metric (6) is written in Fefferman- Graham coordinate z. It is useful to express it in 
a more canonical form, that will prove usefull for a paralell calculation with the static case. 
Indeed, by introducing the coordinate change 

the metric (6) is now 

„2t^2d2 p2 J 2 

ds 2 = j [-f(u)dr 2 + r 2 dy 2 + dx 2 ± )} + 

R 2 R 2 r _1 

+ — T" 2 dr 2 ; drdu , (9) 

9 3 u^ftuj 



3 



with / = 1 — u 2 . We may ingnore the last two terms since the perfect fluid geometry is valid 
only in the scaling limit r — > oo and v,u — > constant. A further transformation through, 



yields 



ds2 = ^ 



u 



t/to = ^(r/r )^ , 



-f(u)dt 2 + ^t 2 dy 2 + ^jdxl 



(10) 



+ 



R 2 du 2 
4/(u) ii 2 



(11) 



This form is similar to the canonical black hole metic except for the time dependence of the 
spacial parts. Also, t and u are not the same variables as in the static case. They are related 
to the original variables, x° and z through (7), (8) and (10). In this transformed metric (11), 
the black hole horizon is no longer moving away in the fifth direction but is expanding in 
the y direction and contracting in the transverse direction as time goes on. 

In the background (11), let us consider the small string fluctuations in the transverse 
direction x±, 



5X 1 = Z(t,u) . 



(12) 



The relevant Nambu-Goto action is 



S = 



T V\ 



dt 



du ( — - 

2t 



(ftO 2 4/( M )tt 2 T 2 



U3f(u) 



U2 



(<U) S 



(13) 



where v^A = ^ after subtracting the unperturbed string action. Notice that the action (13) 
is the same as the one in the static black hole metric [20] except for an overall factor of (^) • 
The latter stems from the metric component g x i x i of (13), when we evaluate the induced 
metric in the Nambu-Goto action. The equation of motion for £(t,u) is 



m - fti + 2vr 4 T 7H(l + 3« 2 )cU - 47r 4 T VH 2 ^ = . 
To solve (14) we define a Fourier-like transform 



(14) 



d(t,u) = f^^tHf\ut)^MU^), (15) 
where Hf\uot) is a Hankel function of the second kind, and ^(w) is normalized such that 
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V&uXO) = l. tHf\ujt) is chosen to satisfy the time part of the equation of motion (14) with 
the correct boundary condition (see below). Notice that we have extended the region of t 
from (0,oo) to (—00,00). 

To proceed, we assume the following 'completeness relation' 

-=- / dt tH?\ujt)H?\-uj't) ~ -6(u - to'). (16) 
4 J-00 w 

which will be understood for small u 1 . While this approximation blurs the rigor of the 
arguments to follow, it should nevertheless provide us with an insightful understanding of 
the time scale involved in the relaxation of the diffusion process. These time scales are 
paramount to our understanding of the approach to equilibrium of a strongly coupled quark- 
gluon plasma such as the one at RHIC. 

After separating the time part, the equation for ^^(u) now reads 

3u 2 + 1 tt> 2 

~ ww^W + S7^*" ( ") = ■ (17 > 

where ro = Notice that (17) is of the same form as the one in the static black hole 
metric [20]. Near the horizon the solution behaves as 

^~(1- M ) ±IW / 4 , (18) 

and the minus choice corresponds to the infalling boundary condition. Inserting (15) into 
the action (13) yields the reduced boundary action 

_ 37T 2 v^T 3 t f /(«) 

^boundary - ^ J Ot — ^-£(t, U)O u ^{t, U)\ U=Q 

dcu ~ , A( 37r 2 -s/A7ot \ f(u 



= / ^ 



u 



u=0 



1 The usual completeness of the Hankel transform is in terms of Bessel functions: J^°dt tJ„(u)t)J v (u>'t) = 
—5(u> — oj'), whose origin is the asymptotic form of a Bessel function as an exponential function over \ft. 
There is not true completeness for the Hankel function #(1,2) due to the singularity near zero. The use the 
completeness of the Hankel function is justified for small u> or large times since: l)the dominant integral 
contribution is coming from the large time region; 2)the JP background is justified asymptotically. The use 
of the Hankel function instead of a Bessel function is needed to match the incoming bioundary condition 
below. 
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where we used (16). Following [23] we identify the retarded Green's function, Gr(u), as 



G R {u) = 



37r 2 V\Th, 



m 



u=0 



(20) 



3 Electric Force-Force decorrelation 



Gr(oj) can be calculated analytically only in two limits: uj — > and uj — > oo. In the small 
uj limit, we may expand the solution in terms of to and solve (17) order by order with the 
incoming boundary condition, 



1 c 



ivo 



which gives 



1 + — (- tan" 1 yfa + ln(l + y/uj) 



lim (Trcoth^) Im[^@*_ w (u)9 tt * w (u) 
^-►o V 2 / L a/t/ 



+ 0(m 2 ) , 



(21) 



y/u 



U=0 



(22) 



In the large uj limit, we can use the WKB approximation, which yields 



lim (vrcoth^ hn\£^V^(u)dJ& u (u) 
^-►o V 2/1 \ u 



y/u 



7r tn 



J u=0 



(23) 



in agreement with the zero temperature result [21]. See Appendix A for more details. 

For general u>, we have to resort to numerical methods. We should perform the numerical 
computation and compare it with the analytic result in the small and large uj limits obtained 
above. The strategy is as follows. 2 First we find two independent solutions near the horizon 
(u ~ 1) 



H — 



U},OUt 



(1 - u)- iM ' A 

\ u>,iJ ' 



l-(l-li) 



2tD z 



Si + 4m 



+ 



(24) 
(25) 



Here ^" in is the infalling solution and its complex conjugate is the outgoing solution. Notice 
that these solutions are valid for all to. Near the boundary (u ~ 0), there are two independent 



2 Similar calculations have been done in [25, 26, 21] in other models. Our numerical integration is different. 
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solutions 



oj,0 



= u 



10 



It) 

7 + 280 



M 7 / 2 + 



to 2 ro 4 2 / ro 2 ro 6 > 



(26) 
(27) 



Notice that ^f )0 vanishes, while ^ x goes to unity near the boundary and both solutions 
are real. For the retarded Green's function, we need the wave function near zero satisfying 
infalling boundary condition at the horizon. For this, we take the near-horizon wave-function 
(24) with the correct boundary condition as the initial data and numerically integrate it from 
the horizon to the boundary using (17). The solution is expressed as a linear sum of boundary 
basis and 



(28) 



where A and £ are complex numbers determined numerically. Notice that the right hand 
side goes to A at the boundary while we have to normalize ^ such that it goes to 1 at 
u = 0. Therefore the correctly normalized wave function with correct boundary conditions 
is * w =^l- 1 *J in ( W ): 



which readily yields 



Im 



/(«) 



u 



*-»9A(«) 



= § ImS , 

«=o 2 



rn 

with S = Now the Wightman function G(u) (3) is given by 



G(u) 



(29) 



(30) 



(31) 



To complete the numerical calculation, we note that while A is easily accessible numer- 
ically in (28), "B is not. To resolve this problem, we use the following method. First, by 
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G ^ GCJ.s) 




2 27 



[The long dashed red line: discrete Fourier transform of (a). The short dashed blue line: the 
divergent contribution alone. The solid line: the total result.] 



taking the imaginary part of (29) we get 



ImS 



u 



(32) 



and then we evaluate it at any point, say, u — 1. The only remaining part is the value 
of ^> B (u) at u — 1, for which we need to numerically integrate from the boundary to the 
horizon. We denote the value determined by this procedure by ty^ (u ^>1). Notice A is 
given before by ^^ in {u ^Q>0). Therefore we get the numerical recipe: 



ImS = Im 



\D H . (u = I 



<in(u ^0) ■ M> 



B 

u>,0 



U 



(33) 



where we take e = 10 -6 . 

The result for (31) is plotted in Fig. la, with the large ro asymptotic ^jp- subtracted. 
There is a good agreement asymptotically. The numerical results in time can be obtained 
by using the inverse transformation of (16) 3 



i r°° 

G(t h t 2 ) — — — / duuH^^t^H^i-ut^G^) . 

^ J — oo 



(34) 



3 The same comments footnoted after (16) apply. 
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The time-dependence of the problem excludes time-translation invariance. In terms of the 
relative and CM coordinates 



s = tl -t 2: ^^p, (35) 

the force-force decorrelator is 

G(7, s) « I A_ J ^ e—G{u) , (36) 

v 1 A* 

for t\,t2 1, which is an ordinary Fourier transform in the relative time. For s C J the 
asymptote G{uj) ~ |ro| 3 , yields 



/ 



e -<wa M 3 ~ ~ | ^ — rj- (37) 



We note that this is a good approximation of the Hankel transformation (15) for large 
ti. The result (37) is expected from the conformal dimension which is 4 of the force-force 
decorrelator. Similary for s ~ T 

G(7, s) 1 , (38) 

TYT 2 -s 2 /4 

but it is not reliable due to the nature of our approximation (ti,t 2 ^> 1). In Fig.la, we 
show the regular part of G(uo). In Fig. lb, we show its discrete Fourier transform in dashed 
red versus snT . The short dashed blue line is the divergent contribution alone which is 
dominant at small relative times. The solid line is the total result. The decorrelation time 
follows readily from the long dashed red curve as 

t F ~ \ . (39) 

This result is important as it indicates that from dual AdS/CFT all electric-electric forces 
applied to either heavy or light fundamental probes decorrelate on a short time scale of the 
order of 2/irT in the static but strongly coupled QGP. We also note that the decorrelation 
curve in long dashed-red is stronger than exponential. This time compares favorably with 
the time read from the lowest quasi-normal mode tvf 1 associated to string fluctuations (see 
Appendix B) 

tuf w 2.69 - 2.29i (40) 
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This yields a decorrelation time of order 0.44/T which is comparable to our 0.64/T numer- 
ical estimate 4 

Using the relation (10), 



t F = 5t = (ro/r) 1 / 3 ^ , 



(41) 



this 'static' time translates to a 'moving' time 

(r/r ) 1 /3_ l 



5t 



T(r) 



(42) 



which is the natural time dependent temperature, T(r), in agreement with Bjorken hydro- 
dynamics [10, 6]. 

4 Momentum fluctuation and Diffusion of a heavy quark 

The short force-force decorrelation time assessed above justifies the arguments presented in 
section 2. Again in the time window tp < At < to, the fluctuations in the momentum of a 
massive quark with MD > 1 read 



(Ap(t) 2 ) = ({p(t + At) - p(t)f 



(43) 



For a quark at rest this accounts for its momentum broadening by thermal quicks. Using 
(22), then 



ft+At pt+At pt+ist poo 

(Ap(t) 2 ) = / d*x / dt 2 (F(t 1 )F(t 2 )) « / dT/ dsG(7,s) 

Jt Jt Jt J -co 



t+At poo 



= ln|l + T 



lim fyrcoth^) Im[^^*_ w («)9 tt * w , (u) 
ro->0 \ 2/1 \ u 



u=0 



3 r-„ At 

-7rv / AT ( 3 t 0T , 



(44) 



where we have used the fact that G(7,s) is well localized. From (10), the time dependent 
momentum transfer is 



At 



{Ap(r) 2 ) = 7cV\T^t := «(r)Ar . 



(45) 



4 The same value was reported in [21] for the fluctuation of the trailing string. 
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For Brownian diffusion which is the case here given the short decorrelation time of the electric 
force, this amounts to 



K T 



t/t 



ttVXT 3 (t) . (46) 



with T(r) defined in (42). k(t) is the time-dependent momentum diffusion constant. 

Now consider the diffusion of a nonrelativistic heavy quark in the medium of which 
temperature is cooling down adiabatically. From above it follows that the Langevin equation 
captures the essentials of the equilibration in the diffusion regime. Since the medium is 
expanding, the appropriate description is given by 

= -Vd(t) P (t) + F(t) , (F(r)) = 0. (47) 
where t)(t) is a time-dependent drag coefficient which is related to F(t) by 

" M = fa / WW<*)> = ■ (48) 

One way to derive this is to follow [28], where all arguments are done for small time steps 
At whereby the expanding medium is frozen. 

Multiplying both sides of (47) by x{r) and taking the ensemble average of the product 
one gets an ordinary differential equation, 

d2 : ( x ^ +r] (r)^(x 2 )-2(v(r) 2 ) = 0. (49) 



dr 2 dr 

where the only property of F(t) we need is (x(t)F(t)) = (x(t))(F(t)) = 0, which is a basic 
assumption of the Langevin equation. The time-dependent diffusion rate D{r) reads 

D{r) ee \j-(x 2 ) . (50) 

so that (49) is 

D(t) + t]{t)D{t) - (v(t) 2 ) = . (51) 
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To solve (51) we need two inputs: t](t) and (v(r) 2 ). rj(r) is given by (46) and (48), 



which is similar to the static case [17, 19, 20] except for the time dependent temperature. This 
confirms the adiabatic nature of the expansion in the case of short force-force decorrelations. 
Using the adiabatic form of the equipartition theorem, we have 

< v{rf >= ^ . (53) 
Using (52) and (53), we now have 

D(t) + a t~ 2 ^D(t) - b r- 1 / 3 = 0, (54) 
with a = r/ T^ 3 and b = T t^ 3 /M. The solution is 

D(r) = -r 1 / 3 + J D(0)e- W/3 . 
a 

This result is important as it shows how the diffusion rate for a quark changes in an expanding 
medium. At short times it is D(0) while at large times it asymptotes 

flM = wfe (55) 

which is the result in [20] with an adiabatically changing temperature. In a way, this justifies 
a posteriori the use of the 'completeness relation' (16) for large times 7. The cross over 
between short and long times is exponential and of order 3ar x / 3 = 1 in time, so (55) is 
reached for t/t = (l/3^oTo) 3 - At RHIC r w 1 fm so that t/t l/rj^ in the Bjorken 
phase. Although our arguments rely on large times within the diffusion window (see above) 
the cross over regime can still be approached albeit from above. 



5 Conclusion 

We have analyzed the diffusion of a heavy quark in an expanding and strongly coupled QGP 
using the AdS/CFT construction. Our arguments provide some insights to a truly non- 
equilibrium phenomenon at strong coupling. Our analysis was restricted to an asymptotically 
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Bjorken expanding fluid which is dual to the JP metric. 

In the comoving frame the time-dependent diffusion problem is mapped onto a time- 
independent one, whereby the diffusion is captured in a retarded Green's function with 
proper boundary conditions in the gravity dual space. The Green's function reflects on 
the electric nature of the noise in the rest frame of a massive quark. It is important to 
note that the noise (force-force) decorrelation is short with t F ~ 2/nT for both ligh and 
heavy fundamental quarks. A strongly coupled QGP randomize the electric correlations very 
efficiently. 

If we denote by Tp = MD/T the diffusion time, for times tp < t < tp a diffusion 
regime for massive quarks open up whereby a generic and memoriless Langevin description 
holds. This desciption requires the knowledge of only two underlying moments of the phase 
space distribution: the average shift captured by the drag and the diffusion constant. Our 
construction allows the generalization of these concepts to a time-dependent Langevin de- 
scription pertinent for an expanding fluid. We have found that asymptotic diffusion sets in 
on a time scale r ~ l/Vo- This estimate is reached from the generic character of the Langevin 
description. 

Finally and while we have used the construct in [3] for the retarded Green's function, it 
will be of some interest to check explicitly this using the arguments presented in [30]. 
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A WKB for large ro 

In this appendix we solve (17) for large ro. To use the WKB approximation in [23, 25], we 
need to transform (17) to a Schrodinger type equation. By a change of variable 



VT 



u 



2 



= JL Tr7^*» > ( 56 ) 



u 1 
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(17) transforms to 

d 2 \ r i n t t -5 + 18« 2 + 3« 4 + 4«tt> 2 

d^-^ = ' V ^ im-^ ■ (57) 

where V plays the role of a potential. V will be approximated as follows, 

4 + to 2 



u « 1 
WKB 



16(1 -w) 2 ' 
tt) 2 

4m(1 -w 2 ) 2 ' 



" K0: is?-£- (58) 



and the corresponding solutions are 

1i ~ 1 ; (1-m)2 4 ) (i_ M j 2 +4 ; 



2vMl- 


M 2 ) 


tt> 


2^(1 - 


li 2 ) 


tt) 



m ~ : ^I/I 6 ^ 1 " iy/uto) , ^l^e-^ ro (l + iyfao) . (59) 

Since the solutions are valid only in the limited region we need to tie them to fulfill the 
expected boundary conditions. There are two boundary conditions: One is the incoming 
boundary condition near the horizon and the other is the normalization (0) = 1). The 
incoming solution near the horizon is ~ (1 — u)~^ and it corresponds to = (1 — u)^~ 1J ^ . 
It can be shown that all the first parts of the solutions are connected to each other so that 
the physical solution near u = is 

tf w , = e ^ w (l - iVuto) , (60) 

which gives us 



lim Im 



tt_„(u)0 u tt w , («) - — - . (61) 



u=0 
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B Quasi-normal modes 



In this section we compute the lowest quasi-normal mode following the method presented 
in [29]. By a change of variable 

y=l-u, (62) 
Eq.(17) can be reduced to the Heun equation: 

where y=0(at horizen) is a regular singular point with characteristic exponent {— 
The quasinormal mode is the solution of (63) obeying the incoming boundary condition at 
the horizen y = and the vanishing Dirichlet boundary condition at y — 1. So we choose 
the exponent — % j at y — 0. To match the boundary condition at y — 1, it is convenient to 
transform (63) once more by, 

^(y) = y -T( y -2)-T X (y) . (64) 

Then (63) is reduced to the standard form of the Heun equation. 

<&v) + + ^ + j^) 9 ^ + = ' (65) 

where 

irv 2 tt),., .tt) 

a/3:=_--(l +l ) , 7: =1_2_, (66) 

x 1 , tt) tt) 2 ro 

= — , e=l , q := 2 — . 

2 ' 2 ' y 8 V ; 4 

At |/ = 0, the local series solution corresponding to the zero characteristic exponent and 
normalized to 1 is given by 

oo 

Xo (y) = J2 a ^)y n > ( 67 ) 

n=0 
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where 



a = 1 , ai = ^- , (68) 

a n+2 + A n (to)a n+1 + S n (m)a n = (n > 2) , 
A (m) ■= -i n + !)(2^ + e + 3(n + 7) + q 



B n (to) 



2(n + 2)(n + l + 7 ) 
n 2 + n(7 + 5 + e - 1) + af3 



2(n + 2)(n + l + 7 ) 
At y=l, we get the boundary value using the series (67) 



Xo (l) = $>nM , (69) 

n=0 

To find the quasinormal modes, we need to find the zeroes of (69) in the complex to plane. 
This is done by truncating the series after a large number of terms 

N 

|Xo(l)TH£«»(ttOI 2 = 0. (70) 

n=0 

We search for the minimum of |xo(l) Ar | 2 , and check that the minimum value is zero. 
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